options(width=108)
out=F
# The list of subjects, the order of conditions, and the thresholds are derived from Subjects.xlsx
library(xlsx)
library(ggplot2)
library(plyr)
library(reshape)
library(matrixStats)
library(gridExtra)
library(lme4)
library(lmerTest)
library(BayesFactor)
#library(splines)
db <- '/home/egor/Dropbox/' # on Linux
db <- '/Users/Egor/Dropbox/' # Windows
# db <- '~/Dropbox/' # on Mac
source(paste(db, 'Prog/R/myFunctions/pvalfn.R', sep=''))
# theme for plotting:
alpha <- .7
w <- .56 # proportion width in group plots
themefy <- function(p) {
p <- p + theme_bw() +
theme(panel.grid.minor.x=element_blank(), panel.grid.minor.y=element_blank(),
axis.text=element_text(size=8), axis.title=element_text(size=9),
legend.text=element_text(size=8), legend.title=element_text(size=9),
legend.key = element_blank(), legend.margin=margin(t=-.04, unit='in'),
legend.background = element_rect(fill='transparent'),
plot.title=element_text(face='bold'))
}
cc <- c('#e31a1c','#fdbf6f','#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99')
xLab <- 'Target Peak Time (s)'
yLab <- 'Log Contrast Threshold'
# colLab <- expression(paste('\nTarget\nVelocity (', degree, '/s)', sep=''))
colLab <- expression(paste('Target Eccentricity (', degree, ')', sep=''))
# colLabType <- 'Mask Type'
yLim <- c(-1.5,-0.5)
dodge <- position_dodge(width=0.1)
sumFn <- function(ss, subjStr='subj', xStr='targTpeak', grpStr='targEcc'){
sumSubj <- ddply(ss, c(subjStr, xStr, grpStr), summarise,
mnS=mean(threshMean), se=sd(threshMean)/sqrt(length(threshMean))) #, .drop=F)
sumSubjMn <- ddply(ss, c(subjStr), summarise, mnStot=mean(threshMean)) # total mean across conditions per subj
sumSubj <- merge(sumSubj, sumSubjMn)
sumSubj$normS <- - sumSubj$mnS / sumSubj$mnStot # normalized subject mean
sumSubj$seNorm <- NA
# sumSubj$mnS[is.na(sumSubj$mnS)] <- 0
sumGrp <- ddply(sumSubj, c(xStr, grpStr), summarise,
mn=mean(mnS), se=sd(mnS)/sqrt(length(mnS)),
norm=mean(normS), seNorm=sd(normS)/sqrt(length(normS)))
sumGrp$subj <- 'average'
sumSubj <- rename(sumSubj, c(mnS='mn',normS='norm'))
sumComb <- rbind(sumGrp, subset(sumSubj, select=-mnStot))
sumComb$se[is.na(sumComb$se)] <- 0
sumComb
}
plotAve <- function(pss, subjStr='subj', xStr='targTpeak', grpStr='targEcc',
xlab=xLab, ylab=yLab, collab=colLab, yStr='mn', seStr='se'){
pss$yMin <- pss[,yStr] - pss[,seStr]
pss$yMax <- pss[,yStr] + pss[,seStr]
pss[,grpStr] <- factor(pss[,grpStr])
p <- ggplot(pss, aes_string(x=xStr, y=yStr, colour=grpStr, group=grpStr,
ymin='yMin', ymax='yMax')) +
geom_point(position=dodge, size=1, alpha=alpha) + geom_line(position=dodge, alpha=alpha) +
scale_x_continuous(breaks=c(.5,1,1.5), labels=c('0.5','1','1.5')) +
geom_linerange(position=dodge, show.legend=F, alpha=alpha) +
labs(x=xlab, y=ylab, colour=collab) + #ylim(0,1) +
guides(colour=guide_legend(keyheight=.3, default.unit='inch'))
p <- themefy(p)
}
plotIndiv <- function(pss, subjStr='subj', xStr='targTpeak', grpStr='targEcc',
xlab=xLab, ylab=yLab, collab=colLab, yStr='mn', seStr='se'){
pss$yMin <- pss[,yStr] - pss[,seStr]
pss$yMax <- pss[,yStr] + pss[,seStr]
pss[,grpStr] <- factor(pss[,grpStr])
p <- ggplot(pss, aes_string(x=xStr, y=yStr, colour=grpStr, group=grpStr,
ymin='yMin', ymax='yMax')) +
facet_wrap( ~ subj, ncol=4) +
geom_point(position=dodge, size=1, alpha=alpha) + geom_line(position=dodge, alpha=alpha) +
geom_linerange(position=dodge, show.legend=F, alpha=alpha) +
scale_x_continuous(breaks=c(.5,1,1.5), labels=c('0.5','1','1.5')) +
# scale_y_continuous(breaks=c(0,.25,.5,.75,1), labels=c('0','','0.5','','1'), limits=c(0,1)) +
labs(x=xlab, y=ylab, colour=collab) +
guides(colour=guide_legend(keyheight=.3, default.unit='inch'))
p <- themefy(p)
}
allDataDir <- paste(db,'Projects/mc/data_bv3/mcBv3_xvv',sep='')
dataDirs <- dir(allDataDir)
dataDirs <- dataDirs[grep('xvv',dataDirs)]
colsOfInt <- c('participant', 'dom', 'session', 'mcBv', 'targTpeak', 'targXoff2',
'targV', 'stairStart', 'meanRev6')
df <- data.frame()
dfRevs <- data.frame()
dfIntn <- data.frame()
for(curDir in dataDirs){
print(curDir)
curDf <- read.csv(paste(allDataDir,'/',curDir,'/',curDir,'.csv',sep=''))
curDf <- curDf[,colsOfInt]
curDf$cond <- substr(curDir,19,22)
subjStairs <- dir(paste(allDataDir,'/',curDir,'/',sep=''))
subjStairs <- subjStairs[grep('.tsv',subjStairs)]
for(curStairFN in subjStairs){
curStair <- paste(allDataDir,'/',curDir,'/',curStairFN,sep='')
curRevs <- read.table(curStair, skip=1, nrows=1)
curIntn <- read.table(curStair, skip=4, nrows=2)
curInfo <- readLines(curStair)
curInfoCols <- data.frame(subj=curDf$participant[1], dom=curDf$dom[1],
sess=curDf$session[1], stairStart=curInfo[33],
mcBv=curInfo[41], targTpeak=curInfo[43],
targEcc=curInfo[37], targV=curInfo[31])
curDfRevs <- cbind(curInfoCols, curRevs[,2:9])
nTrials <- ncol(curIntn)-1
curDfIntn <- curInfoCols[rep(seq_len(nrow(curInfoCols)), each=nTrials),]
curDfIntn$trial <- seq(1,nTrials)
curDfIntn$intn <- as.numeric(curIntn[1,2:(nTrials+1)])
curDfIntn$resp <- as.numeric(curIntn[2,2:(nTrials+1)])
rownames(curDfIntn) <- NULL
dfRevs <- rbind(dfRevs, curDfRevs)
dfIntn <- rbind(dfIntn, curDfIntn)
}
df <- rbind(df, curDf)
}
## [1] "mc2_tgT-mcBv3_xvv-cent_p1_s1_2017-10-05_1028"
## [1] "mc2_tgT-mcBv3_xvv-cent_p12_s1_2017-10-10_1604"
## [1] "mc2_tgT-mcBv3_xvv-cent_p13_s1_2017-10-10_1653"
## [1] "mc2_tgT-mcBv3_xvv-cent_p5_s2_2017-10-06_1627"
## [1] "mc2_tgT-mcBv3_xvv-cent_p7_s2_2017-10-06_1432"
## [1] "mc2_tgT-mcBv3_xvv-dyna_p10_s1_2017-10-10_1033"
## [1] "mc2_tgT-mcBv3_xvv-dyna_p11_s2_2017-10-10_1425"
## [1] "mc2_tgT-mcBv3_xvv-dyna_p14_s2_2017-10-13_1137"
## [1] "mc2_tgT-mcBv3_xvv-dyna_p2_s1_2017-10-05_1234"
## [1] "mc2_tgT-mcBv3_xvv-dyna_p6_s1_2017-10-06_1040"
## [1] "mc2_tgT-mcBv3_xvv-dyna_p9_s1_2017-10-06_1549"
## [1] "mc2_tgT-mcBv3_xvv-peri_p13_s2_2017-10-12_1608"
## [1] "mc2_tgT-mcBv3_xvv-peri_p4_s1_2017-10-05_1539"
## [1] "mc2_tgT-mcBv3_xvv-peri_p5_s1_2017-10-05_1640"
## [1] "mc2_tgT-mcBv3_xvv-peri_p7_s1_2017-10-06_1245"
## [1] "mc2_tgT-mcBv3_xvv-stat_p11_s1_2017-10-10_1247"
## [1] "mc2_tgT-mcBv3_xvv-stat_p14_s1_2017-10-12_1240"
## [1] "mc2_tgT-mcBv3_xvv-stat_p2_s2_2017-10-05_1425"
## [1] "mc2_tgT-mcBv3_xvv-stat_p3_s1_2017-10-05_1342"
## [1] "mc2_tgT-mcBv3_xvv-stat_p6_s2_2017-10-12_1035"
## [1] "mc2_tgT-mcBv3_xvv-stat_p8_s1_2017-10-06_1357"
## [1] "mc2_tgT-mcBv3_xvv-stat_p9_s2_2017-10-10_1348"
ds <- rename(df, c(participant='subj', session='sess', meanRev6='thresh',
mcBv='maskV', targXoff2='targEcc'))
ds$targEcc <- round(ds$targEcc / 35,1)
ds$targV <- round(ds$targV / 3.5,1)
ds$maskType <- ''
ds$maskType[ds$maskV==0.01] <- 'static'
ds$maskType[ds$maskV==0.6] <- 'slow'
ds$maskType[ds$maskV==9.6] <- 'fast'
ds$maskV <- round(ds$maskV * 60 / 35, 1)
ds$maskV[ds$maskV<0.05] <- 0
dfIntn$maskV <- round(as.numeric(levels(dfIntn$mcBv))[dfIntn$mcBv] * 60 / 35, 1)
dfIntn$maskV[dfIntn$maskV<0.05] <- 0
ds$condSplit <- ''
ds$condSplit[ds$cond=='stat' | ds$cond=='dyna'] <- 'v'
ds$condSplit[ds$cond=='cent' | ds$cond=='peri'] <- 'ecc'
head(ds)
## subj dom sess maskV targTpeak targEcc targV stairStart thresh cond maskType condSplit
## 1 1 1 1 0.0 0.5 0.8 0.0 -2 -0.9583333 cent static ecc
## 2 1 1 1 1.0 0.5 0.8 4.6 -2 -0.7166667 cent slow ecc
## 3 1 1 1 16.5 1.0 0.8 0.0 -2 -1.6083333 cent fast ecc
## 4 1 1 1 1.0 1.5 0.8 4.6 0 -0.6583333 cent slow ecc
## 5 1 1 1 1.0 1.0 0.8 0.0 -2 -0.8750000 cent slow ecc
## 6 1 1 1 16.5 0.5 0.8 0.0 -2 -1.5916667 cent fast ecc
thresh <- ddply(ds, .(subj,dom,sess,maskV,maskType,targTpeak,targEcc,targV,
cond,condSplit), summarise, threshMean = mean(thresh))
head(thresh)
## subj dom sess maskV maskType targTpeak targEcc targV cond condSplit threshMean
## 1 1 1 1 0 static 0.5 0.8 0.0 cent ecc -0.8916667
## 2 1 1 1 0 static 0.5 0.8 4.6 cent ecc -1.2208333
## 3 1 1 1 0 static 1.0 0.8 0.0 cent ecc -1.1041667
## 4 1 1 1 0 static 1.0 0.8 4.6 cent ecc -1.2333333
## 5 1 1 1 0 static 1.5 0.8 0.0 cent ecc -1.0875000
## 6 1 1 1 0 static 1.5 0.8 4.6 cent ecc -1.2375000
ss <- sumFn(thresh[(thresh$maskType=='static' & thresh$targV==0),])
ssAve <- ss[ss$subj=='average',]
ssIndiv <- ss[ss$subj!='average',]
pTargStatMaskStat <- themefy(plotAve(ssAve)) + ylim(yLim)
(plotIndiv(ssIndiv) + ylim(-2,0))
iss <- dfIntn[dfIntn$maskV==0 & dfIntn$targV==0,]
(p <- ggplot(iss, aes(x=trial, y=intn, colour=targEcc, linetype=stairStart)) +
facet_wrap(subj ~ targTpeak, ncol=6) + geom_point() + geom_line() + ylim(-2,0))
ss <- sumFn(thresh[thresh$maskType=='static' & thresh$targV>0,])
ssAve <- ss[ss$subj=='average',]
ssIndiv <- ss[ss$subj!='average',]
pTargDynaMaskStat <- themefy(plotAve(ssAve)) + ylim(yLim)
(plotIndiv(ssIndiv) + ylim(-2,0))
iss <- dfIntn[dfIntn$maskV==0 & dfIntn$targV==16,]
(p <- ggplot(iss, aes(x=trial, y=intn, colour=targEcc, linetype=stairStart)) +
facet_wrap(subj ~ targTpeak, ncol=6) + geom_point() + geom_line() + ylim(-2,0))
ss <- sumFn(thresh[thresh$maskType=='slow' & thresh$targV==0,])
ssAve <- ss[ss$subj=='average',]
ssIndiv <- ss[ss$subj!='average',]
pTargStatMaskSlow <- themefy(plotAve(ssAve)) + ylim(yLim)
(plotIndiv(ssIndiv) + ylim(-2,0))
iss <- dfIntn[dfIntn$maskV==1 & dfIntn$targV==0,]
(p <- ggplot(iss, aes(x=trial, y=intn, colour=targEcc, linetype=stairStart)) +
facet_wrap(subj ~ targTpeak, ncol=6) + geom_point() + geom_line() + ylim(-2,0))
ss <- sumFn(thresh[thresh$maskType=='slow' & thresh$targV>0,])
ssAve <- ss[ss$subj=='average',]
ssIndiv <- ss[ss$subj!='average',]
pTargDynaMaskSlow <- themefy(plotAve(ssAve)) + ylim(yLim)
(plotIndiv(ssIndiv) + ylim(-2,0))
iss <- dfIntn[dfIntn$maskV==1 & dfIntn$targV==16,]
(p <- ggplot(iss, aes(x=trial, y=intn, colour=targEcc, linetype=stairStart)) +
facet_wrap(subj ~ targTpeak, ncol=6) + geom_point() + geom_line() + ylim(-2,0))
ss <- sumFn(thresh[thresh$maskType=='fast' & thresh$targV==0,])
ssAve <- ss[ss$subj=='average',]
ssIndiv <- ss[ss$subj!='average',]
pTargStatMaskFast <- themefy(plotAve(ssAve)) + ylim(yLim)
(plotIndiv(ssIndiv) + ylim(-2,0))
iss <- dfIntn[dfIntn$maskV==16.5 & dfIntn$targV==0,]
(p <- ggplot(iss, aes(x=trial, y=intn, colour=targEcc, linetype=stairStart)) +
facet_wrap(subj ~ targTpeak, ncol=6) + geom_point() + geom_line() + ylim(-2,0))
ss <- sumFn(thresh[thresh$maskType=='fast' & thresh$targV>0,])
ssAve <- ss[ss$subj=='average',]
ssIndiv <- ss[ss$subj!='average',]
pTargDynaMaskFast <- themefy(plotAve(ssAve)) + ylim(yLim)
(plotIndiv(ssIndiv) + ylim(-2,0))
iss <- dfIntn[dfIntn$maskV==16.5 & dfIntn$targV==16,]
(p <- ggplot(iss, aes(x=trial, y=intn, colour=targEcc, linetype=stairStart)) +
facet_wrap(subj ~ targTpeak, ncol=6) + geom_point() + geom_line() + ylim(-2,0))
grid.arrange(pTargStatMaskStat + #theme(legend.position='none') +
ggtitle('a: static mask, static target'),
pTargDynaMaskStat + ggtitle('d: static mask, dynamic target'),
pTargStatMaskSlow + ggtitle('b: slow mask, static target'),
pTargDynaMaskSlow + ggtitle('e: slow mask, dynamic target'),
pTargStatMaskFast + ggtitle('c: fast mask, static target'),
pTargDynaMaskFast + ggtitle('f: fast mask, dynamic target'),
ncol=2, nrow=3) #widths=c((1-w)/2,w/2))
cent <- function(v){
v <- apply(v,2,function(x){
x <- x - mean(unique(x),na.rm=T)
x <- x / max(x)
})
}
dsc <- ds
centCols <- c('dom','targTpeak','stairStart','targEcc','maskV','targV','sess')
dsc[,centCols] <- cent(dsc[,centCols])
# dsc$maskV <- (dsc$maskV / max(dsc$maskV)) #*2 - 1
dsc$maskV_c <- (dsc$maskV + 1) / 2
dsc$targV_c <- (dsc$targV + 1) / 2
dsc$targEcc_c <- round((dsc$targEcc + 1) / 2,0)
dsc$targTpeak_c <- (dsc$targTpeak + 1) / 2
# dsc$subj <- as.factor(dsc$subj)
# dsc$maskType <- ordered(ds$maskType, levels=c('static','slow','fast'))
dsc$maskType <- factor(dsc$maskType, c('static','slow','fast'))
head(dsc)
## subj dom sess maskV targTpeak targEcc targV stairStart thresh cond maskType condSplit maskV_c
## 1 1 1 -1 -0.546875 -1 -1 -1 -1 -0.9583333 cent static ecc 0.2265625
## 2 1 1 -1 -0.453125 -1 -1 1 -1 -0.7166667 cent slow ecc 0.2734375
## 3 1 1 -1 1.000000 0 -1 -1 -1 -1.6083333 cent fast ecc 1.0000000
## 4 1 1 -1 -0.453125 1 -1 1 1 -0.6583333 cent slow ecc 0.2734375
## 5 1 1 -1 -0.453125 0 -1 -1 -1 -0.8750000 cent slow ecc 0.2734375
## 6 1 1 -1 1.000000 -1 -1 -1 -1 -1.5916667 cent fast ecc 1.0000000
## targV_c targEcc_c targTpeak_c
## 1 0 0 0.0
## 2 1 0 0.0
## 3 0 0 0.5
## 4 1 0 1.0
## 5 0 0 0.5
## 6 0 0 0.0
# ds$targTpeak_c <- ds$targTpeak - 0.5 # [.5, 1, 1.5] -> [0, .5, 1]
# ds$targEcc_c <- cent(ds$targEcc)
pvalfn(lmer(thresh ~ dom + stairStart + sess + maskType * targTpeak_c * targEcc_c +
(1|subj), data=dsc[dsc$targV==-1,]))
## estm se df tVal pVal star
## (Intercept) -1.027 0.063 22.260 -16.315 0.000 ***
## dom 0.005 0.053 10.828 0.092 0.928
## stairStart 0.030 0.009 388.056 3.150 0.002 **
## sess 0.006 0.018 355.291 0.340 0.734
## maskTypeslow 0.184 0.051 388.056 3.646 0.000 ***
## maskTypefast -0.325 0.051 388.056 -6.421 0.000 ***
## targTpeak_c -0.071 0.055 388.056 -1.278 0.202
## targEcc_c 0.068 0.052 389.444 1.308 0.192
## maskTypeslow:targTpeak_c 0.027 0.078 388.056 0.345 0.730
## maskTypefast:targTpeak_c 0.016 0.078 388.056 0.199 0.842
## maskTypeslow:targEcc_c -0.310 0.073 388.056 -4.231 0.000 ***
## maskTypefast:targEcc_c -0.093 0.073 388.056 -1.271 0.205
## targTpeak_c:targEcc_c 0.010 0.080 388.056 0.128 0.899
## maskTypeslow:targTpeak_c:targEcc_c 0.042 0.113 388.056 0.369 0.712
## maskTypefast:targTpeak_c:targEcc_c 0.025 0.113 388.056 0.216 0.829
pvalfn(lmer(thresh ~ dom + stairStart + sess + maskType * targTpeak_c * targEcc_c *
targV_c + (1|subj), data=dsc))
## estm se df tVal pVal star
## (Intercept) -1.075 0.068 20.898 -15.905 0.000 ***
## dom 0.036 0.059 11.908 0.611 0.552
## stairStart 0.029 0.006 752.920 4.465 0.000 ***
## sess -0.027 0.008 758.904 -3.538 0.000 ***
## maskTypeslow 0.184 0.048 752.920 3.839 0.000 ***
## maskTypefast -0.325 0.048 752.920 -6.761 0.000 ***
## targTpeak_c -0.071 0.053 752.920 -1.345 0.179
## targEcc_c 0.062 0.049 753.271 1.262 0.207
## targV_c 0.054 0.049 753.263 1.095 0.274
## maskTypeslow:targTpeak_c 0.027 0.074 752.920 0.364 0.716
## maskTypefast:targTpeak_c 0.016 0.074 752.920 0.210 0.834
## maskTypeslow:targEcc_c -0.310 0.069 752.920 -4.455 0.000 ***
## maskTypefast:targEcc_c -0.093 0.069 752.920 -1.338 0.181
## targTpeak_c:targEcc_c 0.010 0.076 752.920 0.134 0.893
## maskTypeslow:targV_c 0.054 0.069 752.920 0.770 0.442
## maskTypefast:targV_c 0.465 0.069 752.920 6.696 0.000 ***
## targTpeak_c:targV_c 0.035 0.076 752.920 0.458 0.647
## targEcc_c:targV_c 0.048 0.071 752.923 0.677 0.499
## maskTypeslow:targTpeak_c:targEcc_c 0.042 0.108 752.920 0.389 0.698
## maskTypefast:targTpeak_c:targEcc_c 0.025 0.108 752.920 0.228 0.820
## maskTypeslow:targTpeak_c:targV_c -0.018 0.108 752.920 -0.167 0.867
## maskTypefast:targTpeak_c:targV_c 0.056 0.108 752.920 0.516 0.606
## maskTypeslow:targEcc_c:targV_c 0.120 0.101 752.920 1.196 0.232
## maskTypefast:targEcc_c:targV_c -0.009 0.101 752.920 -0.086 0.931
## targTpeak_c:targEcc_c:targV_c -0.065 0.110 752.920 -0.587 0.558
## maskTypeslow:targTpeak_c:targEcc_c:targV_c 0.072 0.156 752.920 0.462 0.644
## maskTypefast:targTpeak_c:targEcc_c:targV_c -0.007 0.156 752.920 -0.042 0.966